/*
 * GenClst_Algorithm.cpp
 *
 *  Created on: May 28, 2011
 *      Author: nzhao
 */

#include "share/GenClst_Algorithm.h"


int * LINK_MAP;

int isLink1(int i, int j);
int	isConnected(int * spinlist, int num);
void Absorb(int* connected, int * nConnected, int* beingAbsorb, int * nBeingAbsorb, int* locallinkmap, int dim);
int	clst_cmp(const void * i,const void *j);

/////////////////////////////////////////////////////////////////////////////////
// member of Single_Cluster
void Single_Cluster::initialization(int n, int * p_idx)
{
	int i;
	nspin = n;
	spinIDX = new int [nspin];
	for(i = 0; i < nspin; i++)
		spinIDX[i] = p_idx[i];
}

void Single_Cluster::set_Subs(int * val, int num)
{
	int i;
	sub_num = num;
	subIDX = new int [num];
	for (i = 0; i < num; i++)
	{
		subIDX[i] = val[i];
	}
}

void Single_Cluster::append_Sup(int n)
{
	int i, * tmp;

	tmp = new int [sup_num + 1];
	for(i = 0; i < sup_num; i ++)
		tmp[i] = supIDX[i];
	tmp[sup_num] = n;

	delete [] supIDX;
	supIDX = tmp;
	sup_num ++;
}

void Single_Cluster::clear()
{
	if(spinIDX != NULL)
		delete [] spinIDX;
	if (subIDX != NULL)
		delete [] subIDX;
}

/////////////////////////////////////////////////////////////////////////////////
// Public member of Cluster_List
void Cluster_List::initialize()
{
	int i;
	clsts = new Single_Cluster [MAX_CLST_NUM];

	bath_spin_num = NBATH;
	nCLST = bath_spin_num;

	start_of_order[0]	= 0;
	end_of_order  [0]	= bath_spin_num;
	num_of_order  [0]	= bath_spin_num;

	int * p_index = new int [1];
	for (i = 0; i < bath_spin_num; i++)
	{
		p_index[0] = i;
		clsts[i].initialization(1, p_index);
		clsts[i].seqIDX = i;
	}

	cout << "Cluster initialization completed. " << endl;
	cout << num_of_order[0] << " clusters of order 1 are generated. " << endl;

	delete [] p_index;
}

void Cluster_List::generate_from_order(int order_n)
{
	int i;
	int st, ed;
	int n_get, current_idx;

	st			= start_of_order[order_n - 1];
	ed			= end_of_order  [order_n - 1];

	start_of_order [order_n]	= ed ;
	current_idx					= ed ;

	num_of_order[order_n]= 0;
	for (i = st; i < ed; i++)
	{
		n_get		=  growth_by( clsts[i], current_idx );

		current_idx += n_get;
		nCLST		+= n_get;
		num_of_order[order_n]+= n_get;
	}
	end_of_order   [order_n] = nCLST;
	cout << num_of_order[order_n] << " clusters of order " << order_n + 1 << " are generated. " << endl;

	qsort(clsts+start_of_order [order_n], num_of_order[order_n], sizeof(Single_Cluster), clst_cmp );

	for (i = start_of_order [order_n]; i < end_of_order   [order_n]; i++)
	{
		clsts[i].seqIDX = i;
	}
}

void Cluster_List::find_sub_clusters()
{
	int i, j, n,res_n;
	for (i = 0; i < CALCULATED_ORDER; i++)
	{
		cout << "searching subsets of order " << i+1 << "... ..." ;
		for (j = start_of_order[i]; j < end_of_order[i]; j++)
		{
			n		= clsts[j].my_nspin();
			res_n	= find_Sgl_subsets(clsts[j], n);
			clsts[j].set_Subs(temp_search_res , res_n);
		}
		cout << " completed! " << endl;
	}
}

void Cluster_List::clear()
{
	int i;
	for (i = 0; i < nCLST; i++)
	{
		clsts[i].clear();
	}
}

/////////////////////////////////////////////////////////////////////////////////
// Private member of Cluster_List
int Cluster_List::growth_by(Single_Cluster mom, int current_n)
{
	int i;

	int son_idx, m_spin, seq_id;
	int * p_idx;

	son_idx = 0;
	for (i = 0; i < bath_spin_num; i++)
	{
		p_idx = legal_and_accept_element(mom, i);
		if( p_idx != NULL )
		{
			m_spin	 = mom.my_nspin() + 1;
			seq_id	 = current_n + son_idx;

			clsts[ seq_id ].initialization(m_spin, p_idx);

			son_idx ++;
		}
	}
	return (son_idx);
}

int * Cluster_List::legal_and_accept_element(Single_Cluster mom, int idx)
{
	int i, j;
	int mom_dim,
		num_less, num_greater,
		is_legal = 0;

	int * mom_idx;
	mom_dim  = mom.my_nspin();
	mom_idx = mom.mySpinIdx();

	num_greater = 0;


	int * p_res;
	p_res = NULL;


	if ( idx >= mom_idx[0] )
		return p_res;

	for ( i = 0; i < mom_dim; i++)
	{
		if ( idx == mom_idx[i] )
			return p_res;
	}
	for ( i = 0; i < mom_dim; i++)
	{
		if ( isLink1(mom_idx[i], idx) )
		{
			is_legal = 1;
			break;
		}
	}

	if( ! is_legal )
		return p_res;

	//////////////////////////////////////////////////////////////////////////
	// pass the legal test

	for (i = 0; i< mom_dim; i++)
	{
		if (mom_idx[i] > idx)
			num_greater ++;
	}
	num_less = mom_dim - num_greater;

	for (i = 0; i < num_greater; i++)
	{
		temp_idxlist[i] = mom_idx[i];
	}
	temp_idxlist[num_greater] = idx;

	for (i = 0; i < num_less; i++)
	{
		for (j = 0; j < num_less; j++)
		{
			if (j < i)
			{
				temp_idxlist[num_greater + 1 + j] = mom_idx[num_greater+j];
			}
			else if (j > i)
			{
				temp_idxlist[num_greater + 1 + j - 1] = mom_idx[num_greater+j];
			}
		}
		if (isConnected(temp_idxlist, mom_dim) )
//                if (isConnected(temp_idxlist, mom_dim) )
		{
			return p_res;
		}
	}

	p_res = new int[mom_dim + 1];
	for (i = 0; i < num_greater; i++)
	{
		p_res[i] = mom_idx[i];
	}
	p_res [num_greater] = idx;
	for (i = 0; i < num_less; i++)
	{
		p_res[num_greater + 1 + i] = mom_idx[num_greater + i];
	}

	return p_res;
}


int Cluster_List::find_Sgl_subsets(Single_Cluster clst, int nspin)
{
	int i, j, res_n,c_j, sub_id;

	gsl_combination * c;

	res_n = 0;
	for (i = 1; i < nspin; i++)
	{
		c = gsl_combination_calloc (nspin,i);
		do
		{
			for (j = 0; j < i; j++)
			{
				c_j = (int) gsl_combination_get(c,j);
				temp_sub[j] = clst.mySpinIdx()[c_j];
			}
			if (isConnected(temp_sub, i))
			{
				temp_clst.initialization(i, temp_sub);
				temp_search =(Single_Cluster *) bsearch( &temp_clst, clsts+start_of_order [i - 1],num_of_order[i - 1],sizeof(Single_Cluster), clst_cmp);
				sub_id = temp_search->seqIDX;
				temp_search_res[res_n] = sub_id;
				assert(res_n < MAX_VOLUME);
				res_n ++;

				temp_clst.clear();

				clsts[sub_id].append_Sup(clst.seqIDX);
			}

		}
		while(gsl_combination_next(c) == GSL_SUCCESS);
		gsl_combination_free (c);
	}

	return res_n;
}

/////////////////////////////////////////////////////////////////////////////////
// Assistant functions

int clst_cmp(const void * clst1, const void * clst2)
{
	int i, nspin, diff = 0;

	nspin = (*(Single_Cluster *) clst1).my_nspin();

	for (i = 0; i <nspin; i++)
	{
		diff = (*(Single_Cluster *) clst1).mySpinIdx()[i] - (*(Single_Cluster *) clst2).mySpinIdx()[i];
		if (diff!=0)
			break;//return(diff);
	}
	return(diff);
}

int isConnected(int * spinlist, int num)
{
    int i, j, res;
    int * local_linkmap =  new int [num*num];

    for(i = 0; i < num; i++)
    {
        for(j = 0; j < num; j++)
        {
            local_linkmap[i * num + j]=isLink1(spinlist[i], spinlist[j]);
        }
    }

    int nConnected = 0,
        nBeingAbsorb  = num;

    int * connected    = new int[num],
        * beingAbsorb  = new int[num];

    for(i = 0; i < num; i++)
    {
        connected   [i] = 0;
        beingAbsorb [i] = 1;
    }

    int nBeforeAbsorb;
    bool isInc = true;
    while(isInc)
    {
        nBeforeAbsorb = nConnected;
        Absorb(connected, &nConnected, beingAbsorb, &nBeingAbsorb, local_linkmap, num);
        isInc = (nConnected > nBeforeAbsorb);
    }
    res = ( nConnected == num ? 1 : 0 ) ;

    delete [] local_linkmap;
    delete [] connected;
    delete [] beingAbsorb;

    return res;
}

void Absorb(int* connected, int * nConnected, int* beingAbsorb, int * nBeingAbsorb, int* locallinkmap, int dim)
{
    int i, j, nc, nb;

    nc = nConnected[0];
    nb = nBeingAbsorb[0];

    for(i = 0; i < dim; i++)
    {
        if(beingAbsorb[i] == 0)
            continue;

        if(nc == 0)
        {
            beingAbsorb[i] = 0;
            connected [i] = 1;
            nc = 1;
            nb--;
        }
        else
        {
            for(j = 0; j < dim; j++ )
                if(connected[j]==1 && locallinkmap[j*dim + i] == 1)
                {
                    beingAbsorb[i] = 0;
                    connected [i] = 1;
                    nc ++;
                    nb --;
                    break;
                }
        }

    }
    nConnected[0] = nc;
    nBeingAbsorb[0] = nb;
}

int isLink1(int i, int j)
{
	int idxij;

	idxij = i >= j ? (i+1)*i/2+j : (j+1)*j/2+i;
	return(LINK_MAP[idxij]);
}

